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ABSTRACT 



Aims. We study the interaction of a low-mass planet with a protoplanetary disk with a realistic treatment of the energy balance by 
doing radiation-hydrodynamical simulations. We look at accretion and migration rates and compare them to isothermal studies. 
Methods. We used a three-dimensional version of the hydrodynamical method RODEO, together with radiative transport in the flux- 
limited diffusion approach. 

Results. The accretion rate, as well as the torque on the planet, depend critically on the ability of the disk to cool efficiently. For 
densities appropriate to 5 AU in the solar nebula, the accretion rate drops by more than an order of magnitude compared to isothermal 
models, while at the same time the torque on the planet is positive, indicating outward migration. It is necessary to lower the density by 
a factor of 2 to recover inward migration and more than 2 orders of magnitude to recover the usual Type I migration. The torque appears 
to be proportional to the radial entropy gradient in the unperturbed disk. These findings are critical for the survival of protoplanets, 
and they should ultimately find their way into population synthesis models. 
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1. Introduction 

Planet formation takes place in protoplan etary disks that are 
comm only found around young stars dBeckwifh & Sargentl 
Il996l) . These disks play a critical role in structuring the forming 
planetary system in terms of masses and orbits, and planet-disk 
interaction is therefore a very important process to study in order 
to understand planet formation. 

Planet-disk interaction affects a forming planet in three 
ways. First of all, the mass that the disk is able to supply to 
the young planet is limited, w hich leads to a maximum gas 
mass that a planet can achieve (jLubow et al.lll999t [Kiev 1999; 
ID'Angelo etaT]|2003bt iBate et alj|2003l) . This limiting mass is 
of the order of a few Jupiter masses (Mj). But also for planets of 
lower mass the disk regulates the mass accretion, which becomes 
especially apparent in t hree-dimensional numerical simulations 
dD'Angelo et al]l2003bl) . 

Second, the disk is able to c hange the orbital radius of th e 
planet through tidal interaction (Gold reich & Tremaind [l980). 
Depending on the mass of the plan et, two modes of migration 
can be distinguished dWardll997l) : type I migration for low- 
mass planets that generate a linear disk response, and type II 
for high-mass, gap-opening planets. Recently a third type of mi- 
gration was put forward, in which strong corotational torques 
force a very fast m ode of migration dMasset & Papaloizoul2003t 
Artymowicz 2004). This Type III migration is very sensitive to 
local density gradients, in terms of migration speed and even 
migration direction, and therefore the structure of the gas disk 
determines the outcome of the migration process. This is unlike 
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Type I and Type II migration, for which all planets move inward 
and only the time scale varies with disk properties. 

Finally, the eccentricity of the planet's orbit may be af- 
fected by the disk. Usually, eccentricity is damped ( Artvm owiczl 
1 19931) . but high-mass planets (M p > 5 Mj) that open deep gaps 
may experience eccen tricity growth ( Papaloizou et al. 2001; 
iKlev & Dirksenl 120061) . However, due to possible saturation 
of corotation resonances planets of M p « Mj may pick 
up significant eccentric ity during the gap formation process 
dSari & Goldreichl2004[). although this require s a slightly eccen- 
tric orbit to start with fcoldreich & Sarill2003l> . 

Observations of the masses and orbits of extrasolar planets 
provide insight in the result of planet-disk interaction. One of 
the most striking results that were obtained is the discovery of 
a new class of planets, the so-called Hot Jupiters. They are gas 
giant planets orbiting very clo se to their central star. It is unl ikely 
that they have formed in situ (Pap aloizou &Terauemll 19991) and 
inward migration through planet-disk tidal interaction has been 
put forward as the explanation for their existence. 

The main problem for this scenario is that the time scale s for 
migration as obtained from analytical (Tanaka et al. 2002j) and 
numeric al dBate et al.ll2003t ID'Angelo et alJl2003bT) arguments 
are usually smaller than or comparable to the total lifetime of 
the disk. This means that theory essentially predicts that all plan- 
ets should fall onto the central star, which poses a problem for 
planet formation theory. If it is already hard to form a single giant 
planet within the lifetime of the disk (Poll ack et al.l ll996). how 
are we going to form a huge amount of them, of which only a 
tiny fraction survives? Clearly a stopping mechanism is required 
that prevents planets from migrating all the way towards the star. 



One possibility is that the planet encounters an inner hole 
in the disk, possibly created by the magnetic field of the 
central star. When there is no disk around to interact with, 
the planet will stop migrating. Another possibility is that the 
planet undergoes Type III outward migration at some point 
dMasset & Papaloizoul 120031) . Finally, direct interaction of the 
disk with a magnetic fi eld may stop th e planet, either by a 
toroidal magnetic field (Terquem 2003) or, for a low-mass 
planet, it s migration may be slowed d own through magnetic tur- 
bulence dNelson & Papaloizoull2004l) . 

However, most numerical hydrodynamical work on planet- 
disk interaction lacks proper treatment of the energy balance , 
with only a few n o table exceptions jD'Angelo et all [2003 a; 
iKlahr &Kleyl 120061: iMorohoshi & Tanakal 120031) . The usually 
adopted locally isothermal equation of state assumes that all ex- 
cess energy generated by compression, viscous dissipation or 
shocks can be radiated away efficiently, thereby keeping the tem- 
perature profile fixed. However, the ability of the disk to cool is 
strongly linked to the opacity, and therefore to the density. 

The effects of temperature structure on planet migra t ion in 
realistic di sks were studied by iJang-Condell & Sasselovl (2003, 
12004 12005b who used a detailed disk model including radiative 
transfer to calculate the torque on the planet in an analytical way. 
They found that temperature perturbations resulting from shad- 
owing and illumination at the surf ace of the disk can decrease 
the migration rate by a factor of 2. iMenou & Goodman! d2004 
also find reduced migration rates in realistic T Tauri a-disks, and 
they point out the importance of sudden changes in the opacity, 
for example near the snow line. 

In this study, we aim at relaxing the isothermal assump- 
tion in hydrodynamical simulations of planet-disk interac- 
tion. The first results on planet migration were presented in 
Paardekooper & M ellemal (12006a), where it was shown that low- 
mass planets may migrate outward in disks with realistic tem- 
perature structures. In this paper, we continue our study of mi- 
gration of low-mass objects, wit h more eye for detail than in 
Paardekooper & Mellema (2006a), and include results on accre- 
tion. 

In contrast with lKlahr & Kleyl d2006l) . we focus on low-mass 
planets that do not open a gap in the disk. For these deeply 
embedded planets temperature effects are the most important. 
However, due to their low mass their gravitational sphere of in- 
fluence is small, and one needs very high resolution locally to 
resolve the atmosphere of the planet. Furthermore, because ra- 
diation is the main cooling agent, we need to perform radiation- 
hydrodynamical simulations. This is a major leap from isother- 
mal simulations, because we are not only adding an equation for 
the gas energy but also the evolution of the radiation field needs 
to be followed. 

The plan of this paper is as follows. We briefly review the rel- 
evant cooling time scale in Sect. 12 effects of radiation in Sect. [3] 
and the condition for convection in Sect.|4] In Sect. [5] we present 
the adopted disk model, and in Sect. [6] we discuss the numerical 
method. In Sect.|7]we present the results, from simple isothermal 
models to the full radiation-hydrodynamical models. We give a 
discussion on the results in Sect. [8] and we conclude in Sect. [9] 



2. Cooling properties 

We take the analytical opacity data from iBell & Linl ([1994), 
which state that for low temperatures (beyond the snow line) the 



Rosseland mean opacity in cm is: 

K = K pT 2 , 



where p is the density in g cm -3 , T denotes the temperature 
(K) and kq — 2.0 10~ 4 is a constant. We discuss the cooling 
properties of the disk as a function of opacity, starting with the 
optically thin regime. 

The time scale for the coupling between gas and radiation is 
given by 



t, 



1 



coupling 



(2) 



where c is the speed of light in vacuum. When £ C ouplmg is larger 
than the dynamical time scale tdyn> which equals the orbital time 
scale of the planet, the gas can not transfer its internal energy to 
the radiation field and therefore the gas will not cool efficiently 
through radiation. Plugging in the opacity of Eq. (JTJ, and assum- 
ing T w 50 K, which is appropriate for the location of Jupiter in 
the solar nebula, we find that the density near the planet should 
satisfy 



p > 10 



-18 g cm 3^ 



(3) 



in order to cool through radiation. At the midplane of a pro- 
toplanetary disk in the planet forming region this condition is 
always satisfied. 

In general, when a region is optically thin and Eq. © holds 
it can cool efficiently because all its thermally emitted photons 
can escape from this region. In a protoplanetary disk the scale 
length of interest is the disk thickness H. For planet-disk inter- 
action specifically, most of the torque on the planet c omes from 
mater ial that is closer than ~ 2H from the planet (Bat e et aT] 
2003). The condition for the disk to be optically thin over one 
pressure scale height is given by: 

kH < 1. (4) 

Putting in our opacity law Eq. ([T|i we can write a condition for 
the density: 

1 



P< 



kqT 2 H 



(5) 



At the location of Jupiter in a typical protoplanetary disk T rj 
50 K and h = H/r = 0.05, where r is the orbital dis- 
tance of Jupiter. Putting in these values, we find that for the 
disk to be optically thin over one pressure scale height p < 
5.1 10~ 13 g cm -3 , while the density at 5 AU in the minimum 
mass solar nebula (MMSN) is given by p = 10~ n g cm~ 3 . 

When a region is optically thick, it can still cool efficiently 
depending on the local conditions. Consider a sphere of radius 
H around the position of the planet. The flux through this sphere 
in the optically thick case can be approximated by 



aT 4 



(6) 



where r = kH is the optical depth over the radius of the sphere 
and a is the Stefan-Boltzmann constant. The internal energy 
density in the sphere is given by: 



P 



p^T 



r-i r-i kqTH(t-iY 



where p denotes gas pressure, Y is the adiabatic exponent, R is 
the universal gas constant and p is the mean molecular weight. 
The cooling time scale is given by the total internal energy con- 
tent of the sphere divided by the total energy that flows through 
the surface of the sphere: 



(1) 
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Putting in the density and temperature appropriate for the loca- 
tion of Jupiter (p = 1CT 11 g cm -3 , T = 50 K and p = 2.4) we 
find that i coo i m 10 idyn, where the dynamical time scale £dyn 
equals the orbital time scale of the planet. Cooling is only effi- 
cient when i coo i < td yn , and therefore we conclude that cooling 
is not efficient for these parameters. 

It is easy to show that the cooling time scale depends on the 
distance to the central star through: 



^cool 



p 2 {r)H 2 {r) 
T(r) ' 



(9) 



When we assume a constant aspect ratio h, which implies that 
the temperature varies with radius as r -1 , and a power law for 
the density with index —3/2 we see that the cooling time does 
not depend on radius. This means that 



'-cool 
^dyn 



10 



AU 



(10) 



This means that at approximately 15 AU in the MMSN t coo \ ss 
tdyn- Inside this radius, cooling will not be efficient. 

Both Eq. (0 and Eq. (0 have a strong temperature depen- 
dence, and therefore it is not possible to determine in advance 
if a disk with an embedded planet that changes the local tem- 
perature and density will be able to cool efficiently. However, 
this analysis shows that for typical densities and temperatures 
including proper cooling mechanisms is essential for planet-disk 
interaction. 

Summarizing, for a fixed temperature of 50 K at 5 AU, we 
can distinguish four cooling regimes: 

- For the lowest densities (see Eq. (0) the gas can not cool 
through radiation. 

- For somewhat higher densities, gas is able to cool through 
radiation but the disk is still optically thin over a distance H. 

- Even higher densities make the disk optically thick over H, 
which means that thermally emitted photons do not immedi- 
ately leave the region that provides the torque on the planet. 

- For the highest densities the disk is optically thick and 



^cool 

ciently. 



^dyii7 which means that the disk can not cool effi- 



Note that only the second regime resembles the isothermal limit, 
while the first regime is not important for low-mass planets that 
remain deeply embedded. In the third and the fourth regime, 
temperature effects may play a large role in determining the total 
torque onto the planet. 



3. Radiative effects 

In general, radiation is not only a cooling agent but it may also 
dynamically change the velocity structure around the planet. The 
effect of radiation on the gas velocity goes through radiation 
pressure, which, in the isotropic case, is given by 



P 



1 



-aT 4 



(11) 



where a is the radiation constant. The gas pressure p is again 
given by the ideal gas law, and we can then write: 
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p a 
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(12) 



For the typical temperature at 5 AU, T — 50 K, we get: 

P _ 1.0 ■ 10~ 17 
P 



P 



(13) 



-li 



radiation 



For a typical midplane density of p — 10~-" g cm" 
pressure is negligible compared to gas pressure. 

This is an important result, because it shows that near a 
deeply embedded planet the dynamical effects of radiation are 
not important. This means that the gas pressure and velocity 
structure will be almost the same as in the locally isothermal 
case. This greatly simplifies the torque analysis: in regions where 
the temperature rises because of gas compression or the release 
of potential energy the density must be lower than in the isother- 
mal case to arrive at the same pressure. Therefore these regions 
will exert a lower torque onto the planet. Any asymmetry in heat- 
ing will therefore show up as a torque asymmetry, which may 
alter the total torque balance on the planet considerably. 

4. Convection 

Heat transport in a planetary envelope can occur through radi- 
ation or through convection. Convection sets in when the tem- 
perature gradient becomes too steep for radiation to transport 
energy towards the surface of the planet. As is well known from 
theory of stellar interiors, this happens when 



dT dT 

dz dz 



(14) 
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where z is the direction towards the surface of the planet and 
the right hand side denotes the adiabatic temperature gradient, 
which can be expressed as: 



dT 

dz 



t r 



a< 1 



h r- 1 



(15) 



We can make a simple estimate of the importance of convective 
energy transport by approximating the actual temperature gradi- 
ent in the planetary envelope as: 



dT 

dz 



T c 



Rr 



(16) 



where T ncD is the temperature of the surrounding nebula, T c is 
the central temperature of the planet and Rr denotes the Roche 
lobe of the planet. The assumption is that the planet connects 
to the disk at a distance given by the Roche lobe, and that 
the disk remains unperturbed from this distance. This is the 
usual assumption in one-dimensional planet formation models 
dPollacket al.lll996h . 

When we further assume that the low-mass planets do not 
change to local pressure scale height in the disk, we can plug 
Eq. (TT3T > and Eq. ( fTol l into Eq. (fl4l to obtain a maximum value 
for the central temperature of the planet: 



T c < 1 



Rr 
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(17) 



For a 5 M® planet embedded in a disk with aspect ratio h = 0.05 
and adiabatic exponent V = 1.4 we find that T c < 2.21 T nc b- 

As will become apparent in Sect. [7] the direction of migra- 
tion critically depends on the radial entropy gradient in the disk, 
or, in other words, whether Eq. ( TBI ) is satisfied in the radial di- 
rection. For a disk that follows a power law in temperature and 



density with indices (3 and a, respectively, the entropy also fol- 
lows a power law with index 



S = f3-{T-l)a. 



(18) 



Our nominal model, with (3 = — 1 and a = —3/2, has a 
negative entropy gradient. Note that because of the rotation 
of the disk, this does not lead to convection. Also, we did 
not find evidence for an in stability due to baroclinic effects 
dKlahr & Bodenheimeril2003l) . 

5. Model design 

Our model consists of three objects: the central star, the planet 
and the disk. Below, we give a short description of each of these 
components. Throughout this paper, we will work in a spherical 
polar coordinate system (r, 9, <fc) that co-rotates with the planet. 

5.1. Central star 

The central star is the main source of gravity and is located in 
the origin of our coordinate system. This makes the system non- 
inertial, the principle upon which the radial velocity searches for 
extrasolar planets are based. Although we take the acceleration 
of the coordinate frame into account, for the low-mass planets 
the effect is negligible. We take the mass of the star to be 1 M . 

5.2. Planet 

Our coordinate frame rotates with the angular velocity of the 
planet, and because we keep the planet on a fixed circular orbit 
it resides at a fixed location on the grid, (r, 6, 4>) — (r p , 7r/2, n). 
The potential of the planet is smoothed over 2 grid cells on the 
highest level of refinement, which is always much smaller than 
the Roche lobe of the planet (see Sect. [SJ. This way, we always 
resolve the potential. The planet is able to accrete matter from 
the disk without changing its dynamical mass. We vary the mass 
of the planet between 0.6 and 300 M ffi for the isothermal models, 
which spans the whole range from the linear regime of Type I mi- 
gration to a gap-opening Jupiter-mass planet. For the radiation- 
hydrodynamical models we focus on a pla net of 5 M m , which is 
well inside the linear regime (see however Masset et al. 2006). 



5.3. Disk 



0.4 r p to r 



The disk is three-dimensional, extending from r 
2.5 r p in the radial direction (r p is the distance from the central 
star to the planet). In the azimuthal direction the computational 
domain is bounded by < <\> < 2tt, while in the polar direction 
we go up to 2.5 pressure scale heights above the midplane of the 
disk: tt/2 - 5H/2r < 9 < tt/2. 

5.3.1. Gas 

The scale height H varies linearly with radius initially, with 
h = H/r = 0.05. The initial temperature profile is found from 
assuming hydrostatic equilibrium in the vertical direction and 
therefore matches the temperature profile of the isothermal sim- 
ulations. The density is also a power law initially, with index 
—3/2. We vary the midplane density at the location of the planet 
to investigate different cooling regimes. Our nominal value is 
10~ n g cm~ 3 , which is appropriate for the location of Jupiter 
in the minimum mass solar nebula. The initial velocities are zero 
in the radial and polar direction, and the Keplerian speed in the 



azimuthal direction, with a correction for the radial pressure gra- 
dient. The equation of state is given by the ideal gas law. 

We model the disk as an ideal, inviscid fluid, and therefore 
its evolution is governed by the Euler equations. These can be 
cast in the following form: 



dW dF r dFg dF d 



dt 



Or 



09 



(19) 



in which W is the state vector, F are the fluxes in the three co- 
ordinate directions and S is the source vector. The state consists 
of the conserved quantities: 



W = r 2 sm9(p, pv r , pvg, pv^,e) T , 



(20) 



where p is the density, v r is the radial velocity, vg and are the 
two angular velocities and e is the total energy density: 

e = \p{vl + r 2 v 2 g + r 2 sin 2 6(1)4, + fl) 2 ) - p$ + ^-j. (21) 

Here is the angular velocity of the coordinate frame, $ denotes 
the potential, p is the pressure and T is the adiabatic exponent. 
The three fluxes are given by: 

F r = r 2 sm9(pv r , pv 2 + p, pv r vg, pv r v<p, pv r w) 1 (22) 
F e = r 2 sin9(pv e , pv r v e , pvg + p/r 2 , pvgv,),, pw e w)^23) 

F,p = r 2 sin9(pv < p, pv^, pvgv^,, pvl + P „ p^w) T (24) 

r z sin 

where w is the enthalpy of the fluid: 



e+p 
w = . 

P 

The source vector S is then: 





prsin 2 6(v<f, + Q) 2 
-2p^vg + sin 9 cos 9p(v<f, - 
-2p^{Q + v<t> ) - 2 cot 9 pvgiVl + v<j> 



(25) 



S = 
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-2p(V rW + Vg-gj + ^gr) - pil^ 



■(26) 



The reas on for not includ ing an anomalous turbulent viscos- 
ity as in iBate et al] (|2003|) and ID' Angelo et al l d2003bl) is that 
the usual a-prescription (IShakura & Sunyaevl II 973) is not a 
good description of the magnetic turbulence that is probably the 
source of viscosity (Balbus & Hawlevlll990l) . The turbulent na- 
tur e of the disk leads to large temporal and spatial variations in 
a dPapaloizou & Nelsonll2003h . Neglecting viscosity altogether 
means that our results are valid in regions of the disk where the 
ionization fraction is not high enough to sustain magnetic turbu- 
lence. 



5.3.2. Radiation 

Planets that do not open gaps (M p <g; Mj) are deeply embedded 
in the disk where the material is very optically thick. It is there- 
fore appropriate to consider radiative transfer in the diffusion 
limit. However, because the density drops more than an order of 
magnitude over a few pressure scale heights above the midplane 
it is necessary to deal with the optically thin regime as well. 

We start by considering a medium at rest, i.e. we set the ve- 
locity of the gas to zero. In the diffusion limit there is essentially 



no angular information in the radiation field, which is then gov- 
erned by a diffusion equation for its energy density E: 



at 3k 



(27) 



where c denotes the speed of light, n is the opacity in cm -1 
and B is the radiative source term. We consider only thermal 
emission from the disk, therefore 



B = aT 



(28) 



where a is the radiation constant and T is the gas temperature. 
In order to conserve total energy an extra source term S Ia d = 
—kc(B — E) appears in the gas energy equation. We work with 
a Rosseland mea n opacity that is typical for the MMSN (see 
iBell & Lin|[T994b . Initially we set E — B in all simulations. 

In the optically thin streaming limit the evolution of the ra- 
diation energy is given by: 



dE 
~dt 



V • (cEh) = kc(B-E), 



(29) 



where n is a unit vector in the direction of propagation of the 
light front. Note first of all that this is an advection equation 
rather than a diffusion equation, and second that there is a strong 
angular dependence of the radiation field. Naively applying Eq. 
( f2Tb in the optically thin case will lead to large errors. 

However, because the most important part of the disk is op- 
tically thick it is still useful to start with Eq. (l27l > but to avoid 
divergences in the optically thin case. Note that avoiding diver- 
gences is different from obtaining the correct solution, which 
would require a method based on discrete ordinates or on a 
Monte Carlo algorithm. Both of these methods have great dif- 
ficulties with the case of a very optically thick disk. 

We write Eqs. (l27l i and (1291 in the more general form: 



dE 

— + V-F R = kc(B - E), 

in which the radiative flux Fr is given by 



cA. 



R 



K 



-VE, 



(31) 



where A is called the flux limiter. The definition of F r 
in Eq. OH in terms of E removes the need for solving 
a differential equation for Fr. Closing the radiative mo- 
ment equations this way is calle d flux-limited diffusion (FLD) 
(iLevermore & Pomraning|[l98ll) . The flux-limiter A deals with 
the transition from the diffusion limit to the streaming limit. The 
exact functional form of A also implicitly determines the angular 
dependen ce of t he r adiation field. For an overview of different 
forms see We have used 



A = <^ 3+V9+I(rp 



10 R+9+V1&0 R+Sl 



for < R < 2; 
for R > 2, 



where 

R= l\™\ 
k E 



(32) 



(33) 



but the results do not depend sensitively on the adopted form of 
A, mainly because of the high optical depth of the disk near the 
planet. 

Flux-limited diffusion typically give s poor results in the 
streaming limit dHayes & NormarJl2003h . but for our case of 



a deeply embedded planet it is a reliable and relatively cheap 
method for doing radiative transfer. 

When the velocity of the gas is taken into account Eq. (f30b 
needs to be extended in two ways to account for advection of 
radiative energy and momentum transfer between gas and radi- 
ation. The full equation for the evolution of the radiative energy 
then reads: 



dE 
~dt 



V • (F R + v E + v ■ V) = kc(B —E) v ■ F R , (34) 

c 



where V is the radiative stress tensor, given by 
V = ET Edd . 

The Eddington tensor Tgdd is calculated using 

^Edd = ^ [(1 - /Edd)I + (3 ./Edd - l)nn) , 



(35) 



(36) 



where n = V E j '|V E\. The Eddington factor /Edd can be cal- 
culated from the streaming factor / s = |.Fr|/c-E: 



h 



Edd 



1/2 



1 6 " s 



for < f B < 0.4; 



/ s ) 2 + / 2 for 0.4 </„<!. 



(37) 



The extra source vector in the gas flow equations S ra d is 
given by: 



\ 



*rad 



o 

K Tp 

\-kc(B-E)J 



The total source term for the gas is given by S + S r 



(38) 



ad • 



(30) 6. Numerical method 



We use the technique of operator splitting to evolve the hydro- 
dynamic part and the radiative part separately every time step. 

6.1. Hydrodynamics 

We numerically solve the flow equations for the gas us- 
ing a three-dimensional versio n of the RODEO method 
(Paardekooper & Mellema 2006b). For isothermal simulations 
the necessary eigenvalues, eigenvectors and projection coeffi- 
cients ar e simple generalizations of the two -dimensional analogs 
given in Paardekoo per & Mellemal (l2006bl) . For simulations in- 
cluding an energy equation, the method is equivalent to the 
non-relativistic version of t he general relativistic method of 
lEulderink & Mellemal (1 1995I) . 

In short, RODEO evolves the Euler equatio ns using an ap- 
proximate Riemann solver (the Roe solver, see lRodll98ll) . to- 
gether with stationary extrapolation to account for the geomet- 
rical and hydrodynamical source terms (Eq. (|26*]|). The Coriolis 
forces are treated in an e xact fashion to enforce exact angular 
momentum conservation (lKleylll998L lPaardekooper& Mellema 
l2006bl) . The extra source terms due to interaction with the radi- 
ation field (Eq. (l38l ) are treated as external source terms and are 
integrated separately (see below). 

A module for adaptive mesh refinement (AMR) is used to 
obtain high resolution close to the planet. Because the planet 
remains at a fixed location on the grid, this treatment amounts to 
having a nested grid. 



Accretion onto the planet is handled in the same way as in 
iD'Angelo et all ((2003b): every time step At the density near the 
planet is reduced by a factor 1 — fAt. The area from which 
we take away mass has a size r acc times the Hill radius of the 
planet. For all accretion r esults in this paper, we have used / = 
5/3 an d r acc = 0.1 (see ID'Angelo et alj|2002h . iKlahr & Klevl 
(2006) used, in addition to runs with this accretion prescription, 
a procedure that conserves total energy locally by adding the 
kinetic and potential energy that is accreted onto the planet to the 
internal energy of the gas. Because we solve for the total energy 
of the gas directly, our accretion prescription does conserve total 
energy. This means that in the accretion region the gas is heated 
by the accretion process. 

When calculating the torque on the planet, we exclude ma- 
terial that re sides within the H ill sphere of the planet. This was 
also done in lBate et al.l (12003). Although material inside the Hill 
sphere may exer t a torque on the planet , which may lead to mi- 
gration reversal (D' Angelo et al. 2003b), this torque is very sen- 
sitive to numerical resolution and the details of the accretion pro- 
cess. Moreover, this is also the region where self-gravity should 
play an important role, which is not included in our models. 
Therefore we focus on torques due to material outside the Hill 
sphere of the planet. We anticipate that for the low-mass planets 
that we consider in this paper this will not affect our results. 

6.2. Flux-limited diffusion 

Because of the intrinsic short time scale for radiative effects, 
which is due to the enormously fast communication speed c 
compared to the sound speed, it is necessary to integrate the 
equation for E in an implicit way. 

First of all, cons ider the gas coupling term . We follow the ap- 
proach described in [Haves & Norman] ( 20031) . The gas is heated 
by the radiation field according to 

de 

— = kc(E - B), (39) 

where e = p/(T — 1) is the internal energy of the gas. Upon 
differencing both sides we obtain: 



e n+1 - e™ = AtKc(E n+1 - B 



(40) 



where X n denotes quantity X evaluated at time n. We approxi- 
mate B n+1 using a Taylor series expansion: 
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and we spell out the internal energy as 

e n+l = P n+1 P nTn+1 



(41) 



(42) 
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n+l 



r-i r-i 

in which 1Z is the universal gas constant and T is the adiabatic 
exponent. 

Putting Eq. d42l and (f4TT > into Eq. d40b we can solve for 

rpn+1 . 

e" + AtKcE n+1 + ZAtncB n 
e n + AAtKcB n 

This expression is used in Eq. fiTt to evaluate B n+1 . We 
are now in a position to difference Eq. (l34l i (for details see 
lHaves & Normanl l2003), which gives us an algebraic equation 
for E™^~1 at grid position (i, j, k): 

di,j,k+lE" jt k+l + d i-l,j,k E i-lJ,k + d i,j-hk E i,j-l,k' 
A 77>ri+ 1 r\n 



(44) 



for certain dij^k- Note that we use only a seven point finite dif- 
ference stencil. This is cheaper than taking into account all 26 
neighboring cells, while we found no noticeable difference in 
the result. 

We can write Eq. d44l as a linear system equation: 



VE 



n+l 



= Q n , 



(45) 



in which V is a N x N matrix and E n+1 and Q n are vectors 
with N elements. Here, N is the total number of grid points. 
Fortunately, this huge matrix T> is sparse; it has only seven non- 
zero entries every row. For such matrices efficient iterative linear 
sol vers have been deve loped in the literature. For an overview 
see lBarrett et al.l(ll994l) . 

The most basic characteristic of a sparse linear solver is 
whether it is a stationary method or not. Stationary methods for 
solving Eq. fi31 l can be cast in the form: 



E 



n+l 



AE"+_\ 



C, 



(46) 



where E™^ 1 is the approximation to E n+1 after m iterations. 
The matrix A nor the vector C depends on m for a station- 
ary method. Examples are the Jacobi method, the Gauss-Seidel 
method and su ccessive over-relaxa tion (SOR), which was used 
for example in IKlahr & Klevl fc006). The main disadvantage in 
using SOR lies in the fact that the convergence behavior of the 
method depends on a relaxation parameter uj, for which the op- 
timal value is problem-dependent. 

Non-stationary methods are a more recent development, and 
they differ from stationary methods in that the structure of the 
update step changes with each iteration. A simple case is given 
by Eq. (|46| | with a matrix A and vector C that do depend on 
in. Examples of non-stationary methods are the conjugate gra- 
dient method (CG), the generalized minimal residual method 
(GMRES), the quasi minimal residual method (QMR) and the 
bi-conjugate gradient stabilized method (Bi-CGStab). 

The optimal method for a given problem depends on the 
structure of the matrix T>. CG, for example, is only applicable to 
symmetric positive definite matrices, while we know in advance 
that due to opacity gradients the matrix T> will not be symmetric. 
One could try CG on the matrix V T T>, where T> T is the transpose 
of T>, but this is not a good idea because the condition number of 
T> T T> is the square of the condition number of T>, which would 
leads to very slow convergence of the method. 

Another difference amongst the various methods is the mem- 
ory requirement. When dealing with large matrices this can be 
a real bottleneck. For GMRES it can be shown that the optimal 
solution is reached within a finite number of iterations, but these 
iterations become more and more expensive because the method 
needs information from all previous steps. To avoid the enor- 
mous storage requirements one could decide to restart GMRES 
after a fixed number of iterations, but this introduces a free pa- 
rameter in the method that needs to be tuned to the problem. 
Finally, QMR is more expensive in memory than Bi-CGStab 
in terms of memory because it needs more auxiliary variables. 
Therefore we have chosen Bi-CGStab as our method. 

One potential disadvantage of Bi-CGStab is that it needs four 
inner products every iteration. An inner product is a global quan- 
tity, and therefore it requires communication between different 
processors, in addition to the boundary conditions. They act as 
synchronization points, and therefore load-balancing is critical 
in this part of the numerical method. We have implemented a 
dynamical load-balancing module, in which during execution of 
the code the processors exchange part of the computational do- 
main to optimize the domain decomposition. Note that due to 
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Fig. 1. Strong scaling test for the three-dimensional radiation- 
hydrodynamical code, shown as the total execution time (rela- 
tive to the execution time using four processors) as a function of 
the total number of processors. The solid line indicates perfect 
scaling, i.e. running on twice as many processors results in half 
the original execution time. Filled circles indicate the hydrody- 
namics part of the code; open triangles the radiation part. Note 
that both parts of the code scale equally well up to 64 processors. 



the AMR it is not true that every processor should get an equal 
part of the computational domain, because at an AMR boundary 
more work is needed (i.e. interpolation of data, flux correction). 

We have preformed a strong scaling test on a low-resolution 
disk with a low-mass planet. The main grid consisted of 128 ra- 
dial cells, 384 azimuthal cells and 8 meridional cells, on top of 
which we put 2 levels of AMR. The simulation was run on dif- 
ferent numbers of processors (n proc ) for a total of 5 orbits of 
the planet, which corresponds to approximately 7500 time steps. 
The resulting total execution times, relative to the case of 4 pro- 
cessors, is shown in Fig.Q] Focusing on the hydrodynamical part 
first (filled circles), we see that the code scales perfectly up to 
Hproc = 8. Due to memory requirements the simulation could 
not run on less than 4 processors, but the perfect scaling from 
4 to 8 processors indicates that the code should scale very well 
from 1 to 4 processors as well. 

When we increase n proc beyond 8, the scaling is still very 
good, although not perfect. This is always the case at some point 
in a strong scaling test, because the amount of data that needs 
calculation stays the same, while the amount of data that has to 
be communicated increases with the number of processors. At 
some point increasing n proc does not result in a reduction in ex- 
ecution time anymore. Fortunately, for up to 64 processors we 
still measure a significant speed-up. This gives us confidence 
that for the high-resolution runs (that have 8 times more com- 
putational cells in the main grid) the scaling will be good up to 
8 x 64 = 512 processors. The simulations presented in this pa- 
per were run on 126 processors of an SGI Altix 3700 CC-NUMA 
system with Intel Itanium-2 processors of 1.3 GHz. No use was 
made of the shared memory capacities of this system. 

An interesting point in Fig. Q] is that the radiation module 
scales equally well as the hydrodynamics, even though the nu- 



merical methods for both parts of the physical problem differ 
significantly. This is an indication that for this problem only a 
small number of iterations is needed in the linear solver to reach 
convergence. 

6.3. Boundary conditions 

The planet excites waves in the disk, which subsequently prop- 
agate to the boundaries of the computational domain. To mini- 
mize reflection of outgoing waves we employ the non-reflective 
bounda ry conditions as outlined in iPaardekooper & Mellemal 
(2006b) for the hydrodynamic variables. At the midplane we 
use reflective boundary conditions. No reflected waves were ob- 
served near the radial and upper meridional (6 = tt/2 — hh/2) 
boundaries. 

For the boundary conditions for the radiation energy density 
we take a different approach, because the initial condition for 
the gas temperature does not correspond to a stationary solution 
for E. In fact, due to the different opacity regimes at different 
heights above the midplane of the disk it is difficult, if not im- 
possible, to obtain a stationary solution for E. In the case where 
the disk is optically thin throughout, the radiation energy density 
should vary as E oc r~ 2 , which means that matter in equilibrium 
with this radiation has a temperature that varies as T oc r -1 / 2 . 
In vertical hydrostatic equilibrium the relative scale height h is 
then proportional to r 1 / 4 , while the initial condition for the gas 
is a constant h. However, all disk models quickly (within a few 
orbits of the planet) adapt to a new equilibrium corresponding to 
the temperature structure set by the radiation field. 

Real protoplanetary disks are heated internally by viscous 
dissipation (in regions where magnetic turbulence operates) and 
externally by radiation from the central star. Irradiation ef- 
fects may play a l arge role in structuring protoplanetary disks 
dDullemondll2000t iDullemond et alj|200lh . but a detailed treat- 
ment of these effects is beyond the scope of this paper. However, 
in absence of heatin g sources the disk cools down rapidly (see 
Klahr & Kley 2006), which will change the torque on embedded 
planets. 

To account for irradiation in an approximate way we set 
the radiation energy density at the inner boundary to a constant 
value, which, for an optically this disk, would correspond to 
h = 0.05 at the location of the planet. This way, we can com- 
pare the radiation-hydrodynamical simulations directly to their 
isothermal counterparts. 

At the midplane we use reflective boundary conditions, while 
at the outer radial boundary we set the ghost cells such that 
Fn.r r 2 , which corresponds to an outflow boundary. At the 
6 = tt/2 — 5h/2 boundary we apply a similar procedure. This 
way, radiation can propagate freely off the grid. 

6.4. Test problems 

The RODEO hydrodynamics module was tested extensively on a 
Cartesian grid with shock tube problems: in all coordinate direc- 
tions separately as well as in two dimensions diagonally. Also a 
genuinely two-dimensional problem, the wind tunnel with step, 
was studied and the results showed that the internal Roe solver is 
able to deal with the complex shock s tructures that arise in this 
problem (see also lMellema et al]ll99lh . On the two-dimensional 
disk-planet problem, RO DEO was compared to other methods in 
Ide Val-Borro et all (120061) . 

The radiation module was tested with one- and two- 
dimensional diffusion problems with a simplified opacity law 
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Fig. 2. Density and velocity fields around an accreting 5 M e planet after 10 orbits. Top panels: slice through z = (9 = tt/2), 
middle panels: slice through y = {(j) = ir), bottom panels: slice through x = — 1 (r = 1). Inside the white contour the velocities 
are supersonic. 



Table 1. Overview of model parameters. From left to right, the columns indicate: planet mass in M ffl , value of the adiabatic exponent 
(T = 1 means a locally isothermal equation of state), viscosity parameter a, radiation flag, accretion flag and the number of 
refinement levels. When n\ cv levels of refinement are used, the local resolution within a distance of 2H from the planet is a factor 
2«iev higher than the base resolution of Ar = 0.0082 r p . For the local models, the indicated resolution is valid for the whole 
computational domain. 
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(k = constant). See IStone et ail (I 19921) for a description of 
the diffusion tests. The numerical results were indistinguishable 
from the analytical solution. 

Another test for the radiation solver is an optically thin cir- 
cumstellar disk. Although FLD is in principle not very well 
suited for optically thin problems, it should nevertheless let E 
diffuse towards the correct solution E oc r~ 2 . Starting from the 
initial condition E oc r~ 4 (for h = constant, Tar -1 , therefore 
E = aT 4 oc r~ 4 ) the numerical solution agrees with E oc r~ 2 
within one orbit of the planet after which it remains stable. This 
is an important test, because it shows that we can compare the 
radiation-hydrodynamical results in the optically thin limit di- 
rectly to the isothermal runs. 

Finally, we considered an optically thick disk with a constant 
opacity. In this case, the equilibrium solution is given by E oc 
r _1 . Again, starting from the initial condition for constant h the 
numerical solution agrees with E oc within a few orbits of 
the planet, depending on the exact value of the opacity. 

7. Results 

In the previous sections we have described the physical 
and numerical setup for the three-dimensional radiation- 
hydrodynamical simulations. However, we now move from the 
simple, 2D, viscous simulations toward the full models in two 
steps. The first step consists of 3D , is othermal simulations tha t 
are compared to lBate et al.1 (l2003h and lD'Angelo etalJ (l2003bh . 
After that, we will include the energy equation but treat the 
cooling in a very simplistic way by setting the adiabatic expo- 
nent r close to unity. Finally, we will present the full radiation- 




Fig. 3. Isovolume for a Mach number of 1 . Inside the solid struc- 
tures the velocity relative to the planet is supersonic. 



hydrodynamical results. In Table Q] we give an overview of the 
simulations discussed in this paper. 




Fig. 4. Density and velocity fields around an accreting 5 M ffi 
planet after 10 orbits. The slice is taken at a distance of 0.005 r p 
from the mid plane of the disk. 



Fig. 5. Slice through = 7r for a 0.6 M ffi planet after 10 or- 
bits, showing density and velocity. For this low-mass planet all 
velocities are subsonic. 



7.1. Isothermal models 

The flow structure around planets in three-di mensional i sother- 
mal disks has been described before dBate et alJ 2003; 
iD'Angelo et ai1 l2"003b). Therefore we limit the discussion to a 
few interesting remarks on the 3D-structure of the accretion flow. 
The disk model we use differs from the one described in Sect. [5] 
in two ways: we use a locally isothermal equation of stat e, and 
a kinematic a-type viscosity dShakura & Sunvaevlll973l) with 
a = 0.004 at t he location of the planet . Essentially this model is 
the same as in ID'Angelo etalJ d2003bl) . 

In Fig. [2] we show density slices along the three coordinate 
axis for a 5 M ffl planet after 10 or bital periods. B y this time the 
flow has reached a steady state (Bat e et alJl2003h . For a low- 
mass planet like this the spiral wave features in the density are 
not strong, but they appear clearly in the velocity field. 

From the middle and bottom panels of Fig. [2] we see that 
the material accreting onto the planet originates predominantly 
above and below the planet. As the material enters the Hill 
sphere of the planet, which for this mass is 0.017 r p large, it 
is accelerated to supersoni c velocities. This was also observed 
by |d5 Angelo et al. (2003b), albeit not for a planetary mass as 
low as 5 M^g. Interesting is the equatorial outflow of mass that 
is apparent in Fig. |2] As the accreting material rains down on 
the planet it squeezes the envelope which forces some material 
to be expelled from the envelope in the equatorial plane of the 
disk. Even before leaving the Hill sphere of the planet this ma- 
terial again reaches supersonic velocities. It is important to note 
that no vertical hydrostatic equili brium is established ne ar the 
planet. This was also observed by lD' Angelo et al.l d2003 b). even 
for non-accreting planets. The appearance of outflows does not 
depend on the boundary at 9 = ir/2: we found the same behavior 
for a disk extending from 9 = tt/2 - bh/2 to 9 = n/2 + bh/2. 
Therefore, this outflow is not linked to the boundary conditions. 

In Fig. [3] we show the three-dimensional structure of the su- 
personic accretion and subsequent equatorial outflow. The solid 



volume represents the region where the Mach number of the 
flow, relative to the planet, is larger than one. For large radial 
distances to the planet the flow is always supersonic due to the 
Keplerian shear; this leads to the two large pizza boxes on both 
radial sides of Fig. [3] The supersonic accretion velocities ap- 
pear as the two small spheres just above and below the planet. 
The equatorial outflow is quickly turned into spiral wave features 
by the Keplerian shear. Note that the outflow is confined to the 
midplane of the disk. This is illustrated further in Fig. |4] where 
we show a density slice parallel to the midplane but at a height 
of 0.005 r p . Note that this is well within the Hill sphere of the 
planet. Even so, the velocity structure differs dramatically from 
that in the midplane. Apart from the unperturbed Keplerian ve- 
locity on both sides of the planet we can also see the horse-shoe 
orbits, as well as a rota ting accretion flow clo se to the planet. 
This was also found by ID' Angelo et all d2003bl) . 

The velocity structure inside the Hill sphere of the planet is 
intrins ically linked to the accretion procedure (ID'Angelo et al.l 
2003b). However, we found that this is not true for the equato- 
rial outflow. It is there always when the vertical accretion flow 
is, albeit with a different magnitude for different accretion pre- 
scriptions. 

For a planet of 0.6 M ffi the flow within the Hill sphere re- 
mains subsonic. The density and velocity field near the planet 
is shown in Fig. [5] Again we can see the vertical accretion flow 
from above and below the planet, some of which is deflected into 
the midplane of the disk where it may leave the Hill sphere of the 
planet. We can identify four circulation patterns near the edge of 
the Hill sphere, which is approximately 0.009 r n large. Some of 
these p atterns can also be seen in the results of ID'Angelo et al.l 
(I2003bh . 

We will now turn to the question whether the differences in 
flow structure we find affect accretion and migration rates of em- 
bedded planets. In Fig. [6] we show the measured accretion rates 
as a function of planetary mass. Note that in this section we do 
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Fig. 6. Accretion rates for two- a nd thre e-dimensional isothermal simulations. Open symbols represent the two-dimensional results 
from Paardekooper & Mellemal d2006bl) . filled symbols denote the three-dimensional results. Left panel: accretion rate in units of 
disk masses per orbit for different planetary masses. Right panel: growth time scales r g = M p /M p in years, where we assume a 
disk mass of 5 Mj. 
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Fig. 7. Migration time scales r m = r p /f p for two- (open sym- 
bols) and three-dimensional (filled circles) isothermal simula- 
tions. As in Fig. [6] we assume a disk mass of 5 Mj. The dot- 
ted line and t he das hed line indicate the analytical results from 
iTanaka et alJ (l2002h for Type I migration in three-dimensional 
and two-dimensional disks, respectively. The horizontal dash- 
dotted line indicates migration on the viscous time scale of the 
disk (Type II migration). 



not limit ourselves to low-ma ss planets in order to co mpare with 
previous numerical results dD'AngeloetaU |2003b). The two- 
dimensional re sults, shown by the open symbols, were discussed 
extensively in Paardeko oper & Mellemal ([2006b) and we show 



them here just for comparison. We have modeled four differ- 
ent planetary masses, spanning a range from deep in the linear 
regime (M p = 0.6 Mg) to well in the non-linear, gap-opening 
regime (M p = 1.0 Mj). 

Comparing the results of two- and three-dimensional simu- 
lations for the two planets with the highest mass, we find very 
good agreement. Especially for the 1 Mj planets this is to be ex- 
pected, because the Roche lobe of this planet is larger than the 
scale height of the disk, which makes the two-dimensional ap- 
proximation applicable. But also for the planet of 0.1 Mj we find 
good a greement, which contradicts the result of lD'Angelo et al.l 
(2003b). However, this discrepancy is small: for this planet, 
ID'Angelo et al.l d2003bl) find a difference of approximately a fac- 
tor of two between 2D and 3D results, while we find a difference 
of 1.3. It is possible that because our method is specifically de- 
signed to handle shocks in a correct way, non-linear effects start 
to play a role for slightly lower planetary masses. Our results on 
migration rates support this suggestion (see below). 

Because the Roche lobes of the low-mass planets are much 
smaller than the disk thickness the two-dimensional approxima- 
tion breaks down. This affects migration as well as accretion, 
but the effect on the accretion rates is the most dramatic. When 
in the two-dimensional approach mass is taken away from the 
location of the planet, what essentially happens is that a whole 
column of gas is accreted onto the planet. In the full three- 
dimensional problem this does not happen, and therefore the ac- 
cretion rates are significantly lower. For the smallest planet we 
consider (M p = 0.6 M®) the difference amounts to more than 
three orders of magnitude (see the right panel of Fig. [6). 

Qualitatively our re sults agree ver y well with the r esults of 
ID'Angelo etaiT(l2003bl) . Interestingly. iBate et al.l (120031) do not 
find an increase in growth time towards the low-mass end, but 
this is probably due to a lack of resolution. Our results show that 
in order to capture significant amounts of gas the mass of a solid 
core needs no be at least several Earth masses. 

Migration rates also differ between 2D and 3D disks 
(ITanaka et al.ll2002T) . Again, the difference is most pronounced 



for low-mass planets, while for planet with masses compared to 
Jupiter the two-dimensional approximation gives valid results. In 
Fig. |7]we show the migration time scales as a function of plan- 
etary mass. For a planet on a circular orbit, the migration time 
scale is related to the torque in the following way: 
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2 in 



(47) 



where L is the angular momentum of the planet and T is the 
torque, which is what we measure during the simulations. 

Again, for the 1 Mj planet the 2D result agrees very well 
with the 3D result, for the same reason as the accretion rates. 
In fact, they approach the migration rate dictated by the viscous 
evolution of the disk (Type II migration). On the l ow-mass end, 
we can nicely reproduce the analytical results from Tanak aet al.l 
(2002) for Type I migration in 2D as well as in 3D. In this 
regime, 3D migration time scales are approximately a factor of 
two larger than for 2D disks. 

In between the two regimes, we find that for a 0.1 Mj planet 
2D and 3D results differ significantly, a lmost an order of magni - 
tude. This is again in contradiction with D'Angel o et alJ d2003b). 
but note that the main difference lie s in the 2D results. These 
differen ces were already discussed in lPaardekooper & Mel lema 
(2006b) and we will not repeat that discussion here. We only 
note that the results for a thre e-dimensional disk agree very well 
with lD'Angelo et all (1200 3b). Also, the departure from the Type 
I analytical result for the 0.1 Mj planet indicates that for this 
planet non-linear effects start to play a role, in accordance with 
the suggestion made above that these effects are also visible in 
the accretion rates. 

Summarizing, we find excellent agreement with analytical 
results of torques in the low-mass regime as well as in the 
high-mass regime. In between, the results are not clear but 
iD'Angelo etafl (120051) showed that this may be a general prob- 
lem for this transition regime. The measured accretio n rates 
agree very well qualitatively with D'Angelo et al.l d2003bh . again 
with the exception of the transition regime between linear and 
non-linear interaction with the disk. 



7.2. Nearly-isothermal models 

As a bridge between the simple, locally isothermal models of 
the previous section and the full radiation-hydrodynamical sim- 
ulations we now consider models that do include the energy 
equation, but in which we treat the cooling in a simplified way. 
Consider the relation between internal energy and pressure for a 
perfect gas: 



p=(T-l)e, 



(48) 



which shows that the adiabatic exponent T is a measure of how 
a change in internal energy affects the gas pressure, and, given 
the density, therefore also the temperature of the gas. For the 
isothermal case, T = 1, pressure and internal energy evolve in- 
dependently and the gas can be compressed indefinitely without 
changing the temperature. To simulate cooling in a very sim- 
plified way one can take a value of T very close to 1, which 
was done by Nelso n~& Benzl (2003a b) for the two-dimens ional 
planet-disk problem. In iPaardekooper & Mellemal d2006bl) we 
commented already that this approach affects planetary accre- 
tion rates to a large extent, because even for T = 1 .001 the planet 
forms a hot bubble around itself that prevents matter from enter- 
ing the Roche lobe. 



In Fig.[8]we show the density and velocity structure around a 
5 M ffi planet after 10 orbits. This figure can be directly compared 
to the isothermal model depicted in Fig. [2] The global properties 
of the velocity field are similar to the isothermal case: inflow 
from above and below the planet and outflow in the midplane of 
the disk. Both inflow and outflow are again supersonic. 

The central density is approximately 2 times lower in the 
nearly-isothermal model compared to the fully isothermal case. 
This is due to the fact that it is not only the density gradient that 
determines the pressure gradient, but also the temperature. This 
means that part of the pressure-support for the planetary enve- 
lope comes from the temperature gradient. We will see below 
that this lower central density substantially affects the accretion 
rate onto the planet. 

Another important difference that can be seen in the top left 
panels of Figs.|2]and[8]is the change in outflow structure. While 
in the isothermal case of Fig. [2] both outflow streams appear to 
be symmetric this is no longer the case when changes in tem- 
perature are taken into account. In the top right panel of Fig. [8] 
we see that the outflow towards positive y appears to be much 
more collimated along the j/-axis than the outflow towards neg- 
ative y, which is responsible for the asymmetry seen in the top 
left panel of Fig. [8] It is important to realize at this point that 
these asymmetries may suffer from resolution effects. Although 
in these simulations we resolve the Roche lobe by approximately 
50 cells, the high-density region inside of the outflow is only re- 
solved by 5 cells. Resolution studies show that the direction of 
the outflow does not depend on resolution, but the magnitude 
of the outflow does. This is because the central temperature and 
density also scale with resolution, as the minimum distance to 
the planet does, and therefore also the depth of the potential well. 
The equatorial outflow therefore critically depends on the con- 
ditions near the center of the planet where the density and the 
temperature reach their highest values. Therefore there exists a 
connection between the conditions deep in the atmosphere of the 
planet and the planet's immediate surroundings. 

In Fig.|9]we show the temperature structure around the same 
planet as in Fig. [H] In this model, the central temperature has 
risen from 50 K initially to more than 100 K. Note that it is 
only a very small region that acquires this higher temperature. 
This reflects the fact that the accreting region is only 10 % of 
the planet's Roche lobe. In this region, the gas is heated by the 
release of potential and kinetic energy of the accreted gas. The 
two streams of outflowing material along the y-axis are clearly 
visible due to their high temperature. Also the collimation of the 
stream in the positive ^-direction is apparent in the top left panel 
of Fig. [9] The outflows remain confined to the midplane of the 
disk until they reach the spiral waves induced by the planet that 
are responsible for upward motion of gas. However, for a low- 
mass planet this upward motion is not strong, and from the top 
right panel of Fig.|9]we see that the hot flow does not reach high 
altitudes. 

Another important effect can be seen in the bottom right 
panel of Fig.|9]that shows the three-dimensional temperature dis- 
tribution around the planet. The solid volume indicates temper- 
atures higher than the initial temperature at the location of the 
planet. The large solid structure at x > — 1 is due to the global 
temperature gradient. To the left of this structure we see the ef- 
fects of the planet, especially the spiral shaped hot outflow. But 
note that the corotational region located at x = —1 is also asym- 
metrically heated: for positive values of y the corotation region 
is locally heated by the planet, while this does not happen for 
negative y. Again, this is a consequence of the asymmetric hot 
outflow that originates deep inside the envelope of the planet. 
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Fig. 8. Same as Fig. 12 but for a nearly isothermal model with T = 1.01. 
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isovolume indicating a temperature higher than the initial temperature at the location of the planet. 



When the pressure remains the same, however, this temperature 
distribution leads to an asymmetric density distribution which 
will affect the total torque on the planet. Because the region be- 
hind the planet is heated compared with the region in front of 
the planet, the density will be relatively low behind the planet 
which leads to a torque that will be positive if the effect is strong 
enough. We will see below that for this model indeed the total 
torque on the planet is positive. 

Before looking at the torques onto the planet we consider the 
measured accretion rates. In Fig.QJJwe show the accretion rate 
as a function of time for three different equations of state: the 
isothermal result discussed in Sect. l7.1l together with two nearly- 
isothermal models. It turns out that the accretion rate depends 
sensitively on the central temperature. This is not surprising, be- 
cause the pressure in near the center of the planet is determined 
mainly by the potential of the planet, which remains the same. 
Therefore, when the central temperature rises by a factor of 2, 
the central density will be a factor of 2 lower which in its turn 



affects the accretion rate onto the planet. From Fig. [9] we saw 
that for 7 = 1.01 the central temperature indeed was a factor 2 
higher than the initial temperature, while from Figs.|2]and[8]we 
see that the central density is lower by approximately a factor 2 
for the model with T = 1.01 compared to the isothermal model. 
In Fig. [10] we see that this change in density is reflected in an 
accretion rate that is lowered again by a factor 2. 

When we take T closer and closer to 1, we expect to retrieve 
the isothermal results. However, as was already mentioned in 
iPaardekooper & Mellema (2006b), even for T extremely close 
to 1 the central temperature of the planet rises significantly. In 
fact, due to the necessary d ivisions by the factor r — 1 (see 
lEulderink & Mellemalll995l) one has to worry about numerical 
stability for F « 1. Therefore the lowest value of T we con- 
sidered is r = 1.005. For this model, the central temperature 
rises to approximately 90 K, which is again reflected in a lower 
accretion rate than for the isothermal model (see Fig.QJJi. 
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Fig. 10. Accretion rates onto a 5 M ffl planet for three different 
equations of state. 




Fig. 11. Total torque onto a 5 M® planet for three different equa- 
tions of state. The torqu es are normalized to the absolute value 
of the analytical result of iTanaka et al.l (120021) . 



In Fig. Q~T]we show the total torque onto the planet for the 
same models as in Fig. [10] The torqu es are normalized t o the 
absolute value of the analytic result of Tan aka et al.l (120021) . and 
as was already shown in Fig.|7]for the isothermal model we find 
very good agreement with the analytical result. This can not be 
said of the nearly-isothermal models, however. Indeed, for both 
values of V we considered we find a positive torque, indicating 
outward migration. But not only is the direction reversed, also 
the magnitude of the torque is surprising: more than six times 
larger than for the isothermal case. The actual value of T has no 
influence on this result. Note that it takes a few orbits longer for 
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Fig. 12. Radial torque profile for a 5 M ffi planet for three differ- 
ent equations of state. The torques are normalized in such a way 
that the total torque for an isothermal equation of state is equal 
to -1. 



the nearly-isothermal models to reach a steady state with respect 
to the torque, which was not the case for the accretion rates. 
This is probably due to the fact that the torque is modified by the 
equatorial outflow. This material has to come close to the planet 
which takes a few orbits. 

It is interesting to find out where exactly this large posi- 
tive contribution to the torque originates. In order to answer this 
question, we show in Fig. [12] the radial torque profile for the 
three differen t equations of state. For the isothermal ca se we see 
the usual (see lBate et alfcOOSH^Angelo et al1l2003bl) contribu- 
tions from the inner and the outer spiral wave, where the inner 
wave exerts a positive torque on the planet and the outer wave a 
negative torque. It is clear that even from inspection by eye that 
the torque due to the outer wave is slightly stronger which leads 
to inward migration. 

The situation is totally different for nearly-isothermal mod- 
els, as can be seen in Fig. [12] For the outer disk the changes 
are marginal, which is not surprising because the structure of the 
outer wave in Figs. [2] and [8] appears very similar. This in con- 
trast with the inner wave, which has a much more open structure 
for the nearly isothermal model compared with the truly isother- 
mal model. This changing wave structure leads to a net positive 
torque on the planet. 

From the temperature plots in Fig. [9] we conclude that the 
outflowing material originates deep within the Roche lobe of 
the planet. This is an indication that accretion and migration 
rates of formi ng planets are c losely linked. This was already 
suggested by D'Angelo et al.l (l2003bl) for the isothermal case. 
However, only when the energy equation is taken into account 
the migration behavior changes this drastically. Although in this 
section we have treated cooling in a simplified way, it provides 
useful insight in what we may expect from the full radiation- 
hydrodynamical models. In the next section, we will consider a 
more realistic treatment of the energy balance by explicitly in- 
cluding radiative cooling. 
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Fig. 14. Same as Fig. [9] but for a non-accreting radiation-hydrodynamical model. 



7.3. Radiation-hydrodynamical models 

The full radiation-hydrodynamical models differ from the 
nearly-isothermal models in two ways: first of all, the adiabatic 
exponent is set to T = 1.4, which reduces the compressibility of 
the gas. Second, the FLD module calculates the energy diffusion 
through radiation. Therefore, the cooling properties of the gas 
depend on the local temperature and density. When cooling is 
not efficient, we expect that due to the higher value of T the ef- 
fects seen in the nearly-isothermal models will be enhanced even 
more. This will certainly happen in the high density envelope of 
the planet, and because we saw in the previous section that tem- 
perature effects deep within the Roche lobe of the planet affect 
the migration behavior we may expect to see similar things in 
the radiation-hydrodynamical models. 

In order to distinguish the effects from the accretion itself 
from the overall radiation-hydrodynamical behavior of the disk, 
we consider two cases in this section: an accreting 5 M® planet 
and a non-accreting planet of the same mass. The non-accreting 
planet is allowed to build a high-density atmosphere because no 
mass is removed from the Roche lobe. However, as was already 



shown bv lD'Angelo eTa l. (2003b) for the isothermal case, even 
in this case the envelope does not reach a state of hydrostatic 
equilibrium in the vertical direction. Material still rains down 
into the Roche lobe from above and below the planet, some of 
which subsequently flows out again. Thus even in the case of a 
non-accreting planet there may exist a close connection between 
the state of the material deep within the planetary envelope and 
the torques due to material outside the Roche lobe. 

In Fig. [13] we show the density and velocity structure for the 
non-accreting planet of 5 M®. The color scale is the same as 
in Figs. 12 and [8] which tells us immediately that the density of 
the planetary envelope is very low compared to the isothermal 
and nearly-isothermal cases. This reflects the effect of the lower 
compressibility of the gas. We will see below that this has dra- 
matic effects on the ability of the planet to accrete gas from the 
disk. 

Looking at the velocity field in Fig. Q~3] we see that in- 
deed there is no hydrostatic equilibrium in the vertical di- 
rection. The same flow pattern appears as in the isothermal 
and nearly-isothermal cases: inflow from above and below and 




outflow along the equator of the planet. Inflow velocities are 
lowe r in this case because the planet does not accrete mat- 
ter dP' Angelo et al.l l2003bl) . The outflow appears much more 
symmetric in all directions compared to the previous nearly- 
isothermal case. Therefore radiative cooling restores the isother- 
mal situation in these regions of lower density than the planetary 
core. However, we will see that still important differences arise 
in the torque on the planet. 

In Fig. [14] we show the temperature structure around the 
planet, on the same scale as in Fig. [9] The most important dif- 
ference with the midplane temperature structure of the nearly- 
isothermal models is the absence of the accretion jets because 
in Fig. [14] we show a non-accreting model. The temperature at 
the center of the planet is also lower; about a factor of 2 higher 
than the initial temperature. This agrees with Eq. ( fT7b for the 
maximum central temperature that can be reached in absence of 
convection. We can also identify the spiral waves in the temper- 
ature structure, because when the gas is compressed in the wave 
the temperature rises. In the dense midplane of the disk the gas 
can not radiate away this energy efficiently, and therefore the 
spiral waves have a higher temperature. 



Inside the equatorial outflow the gas expands, and therefore 
cools down. This effect causes the temperature close to the hot 
core to be lower than the initial temperature. These regions show 
up as dark areas in Fig.[l4j most notably in the top right panel. 

Another difference with the nearly-isothermal models can be 
spotted most clearly in the bottom-left panel of Fig. [14] While 
the heat generated in the envelope of the planet remains confined 
to the midplane of the disk in the nearly-isothermal models (see 
Fig. in the radiation-hydrodynamical models it can diffuse 
to higher altitudes. Because of the large opacity gradient in the 
vertical direction, this is the preferable direction for radiative dif- 
fusion. 

In the bottom right panel of Fig. [T4]we again, just as in Fig. 
[9] show a 3D isovolume of the temperature structure around the 
planet. It is clear that also the spiral wave structure in the tem- 
perature extends further in the vertical direction due to radiative 
heat diffusion. Also it is clear that the strong accretion jets are 
absent. However, what remains is the asymmetric temperature 
structure in the corotation region. There is a bump visible for 
positive y (cj) < </> p , the region behind the planet) while there is a 
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Fig. 16. Accretion rates onto a 5 M 9 planet for radiation- 
hydrodynamical models with two different densities. For com- 
parison, the isothermal accretion rate is also shown. 
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Fig. 17. Same as the right panel of Fig. [6] with the radiation- 
hydrodynamical results (with our nominal initial density of po = 
10 -11 g cm -3 ) included as open triangles. 



small hole for negative y. This asymmetric corotational temper- 
ature structure affects the corotation torque. 

In Fig. [15] we show the temperature and velocity structure 
for the case of an accreting planet. It is immediately clear that 
the temperature is raised considerably in the planetary envelope 
due to the accretion proces. Note that in this case the planetary 
envelope is unstable to convection according to Eq. ( fT7b : the 
heat generated in the accretion process can not be transported 
to the surface of the Hill sphere by radiation alone. Although 
Eq. ( fT7b is of limited applicability because of the assumption 
that the local pressure scale height is the same as for the to- 
tal disk, we can identify heat transport by fluid motion in the 
equatorial plane of the disk. The top left panel of Fig.PT5lshows 
that the outflow velocity pattern is changed with respect to the 
non-accreting model: it seems that heat flows preferentially into 
the outer disk. However, it should be stressed that this observed 
asymmetry may be caused by insufficient resolution at the very 
center of the planetary envelope. Nevertheless, these hot out- 
flows cause significant changes in the temperature and density 
distribution in the immediate surroundings of the planet and are 
therefore potentially important for the migration behavior of the 
planet. The details of the accretion process should therefore be 
subject of further study. 

Also from Fig.[l5]we see that the convective flow is confined 
to the midplane of the disk: vertical velocities in the top right 
and bottom left panel are all directed towards the planet. This 
shows that the vertical thickness of the hot bubble surrounding 
the planet is due to radiative heat diffusion. Its vertical extension 
is roughly the size of the Hill sphere of the planet. 

The formation of the hot bubble has dramatic consequences 
for the measured accretion rate onto the planet. In Fig. [16] we 
show the accretion rate as a function of time for the isother- 
mal case (dotted line) and two radiation-hydrodynamical mod- 
els with two different initial values for the density at the loca- 
tion of the planet. The solid line indicates an initial density of 
po = 1CP 11 g cm~ 3 which is our nominal value. For this density 
we measure an accretion rate that is more than an order of mag- 



nitude lower than for the locally isothermal equation of state. 
The hot bubble surrounding the planet severely limits the gas 
flow towards the planet. Because this bubble is created by the 
accretion process itself, accretion is a self-limited process in this 
case. 

For a lower density the cooling efficiency increases. In fact, 
according to Eq. ( fToT > the cooling time is proportional to the 
square of the density. For a model with a ten times lower den- 
sity (dashed line in Fig.PToTi we measure an accretion rate that is 
almost ten times higher than for our nominal model. This accre- 
tion rate is comparable to the results for nearly-isothermal mod- 
els (see Fig. [Toll. Lowering the density even more did not change 
the accretion rate significantly. This also agrees with the nearly- 
isothermal models, where we saw that even for T as low as 1 .005 
the accretion rate remained far from the isothermal result. This 
can be understood in terms of the high central density required 
to achieve the isothermal accretion rate. From Fig.[2]we see that 
the central density is more than a factor 37 higher than the initial 
density. It is very hard to cool such a high-density envelope ef- 
ficiently, especially when accretion-generated heat is considered 
as well. 

In Fig. [17] we show the growth time scales for the isother- 
mal results (see also Fig. [6J with the radiation-hydrodynamical 
results added as open triangles. The changes with respect to the 
isothermal models are not as dramatic as the changes from mov- 
ing from 2D simulations to 3D simulations, but still it represents 
another order of magnitude slower accretion. Also from Fig. [171 
we see that the effect is slightly less strong for a planet of 0.6 
M® . This is because the potential for this planet is less deep, but 
the decrease in accretion rate is still significant even for this very 
small planet. 

We now turn to the question of planetary migration 
within the radiation hydrodynamical models. From the nearly- 
isothermal models we know that even the migration direction is 
very sensitive to temperature variations, and we expect that be- 
cause for our nominal density cooling is not efficient we should 
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Fig. 18. Total torque onto a 5 M ffl planet for three different radi- 
ation hydrodynamical models. The torq ues are normalized to the 
absolute value of the analytical result of Ta naka et al.l d2002l) . 



recover at least some of the features of the nearly-isothermal 
models. 

In Fig. [18] we show the total torque acting on the planet as 
a function of time, for three different radiation-hydrodynamical 
models. All torques reach a constant value within 10 dynamical 
time scales, even though the cooling time for the high-density 
case is somewhat longer. We checked that the torques did not 
change over several cooling times. Again, as was the case for 
the nearly-isothermal models, the torque needs some more time 
to reach a steady value than in the isothermal case. 

For our nominal, non-accreting model we observe that the 
total torque is again positive, indicating outward migration. Its 
magnitude is comparab l e to t he ordinary Type I torque as cal- 
culated bv lTanaka et alj (120021) . and this is much lower than the 
magnitude of the torque in the accreting, nearly isothermal mod- 
els. It turns out that this is due to the accretion flow. Indeed, 
for our accreting radiation-hydrodynamical model (dashed line 
in Fig. [T8l the magnitude of the torque increases by a factor of 
3 compared to the non-accreting case. It is the asymmetric, hot 
outflow that is responsible for a large part of the positive torque 
in the accreting model. Radiative heat diffusion tends to lower 
the magnitude of the torque, because the flow is spread out more 
evenly around the planet. This effect is not strong, however. 

When we lower the density the magnitude of the torque de- 
creases until at some point it becomes negative. This is because 
the cooling time decreases with decreasing density, and there- 
fore the disk is more and more able to cool towards its initial 
temperature. From Fig. [18] we see that for a 10 times lower den- 
sity the total torque is negative, with a magnitude that is approx- 
imately a fac t or 2 lo wer than the analytical Type I torque from 
Tana ka et all (120021) . As we reduce the density even more the 
torque approaches the analytical result even more, but the scat- 
ter around the equilibrium value increases. This is due to the 
fact that for such densities the disk becomes very optically thin, 
which makes the problem less well suited to treat in the diffusion 
approximation. 



Fig. 19. Radial torque profile for a 5 M ffi planet for three differ- 
ent radiation-hydrodynamical models. The torques are normal- 
ized in such a way that the total torque for an isothermal equation 
of state is equal to —1. 



For accreting models, we do not reach the Type I result, how- 
ever. The results on accretion indicated already that even for very 
low densities there is still a very hot bubble formed around the 
planet, and the subsequent hot outflow comes to dominate the 
torque, which we found to remain positive even for the case of a 
10 times lower density. This is again an indication that accretion 
and migration of embedded planets are closely linked. In order to 
make progress, the detailed process of accretion requires further 
study. 

Nevertheless, even for our non-accreting model we find a 
torque that is positive, which means that the hot accretion flow 
is not the whole story. It is interesting to see where in the disk 
this positive torque arises. In Fig. [19] we show the torque on the 
planet as a function of radius. The low-density case (dash-dotted 
line), giving rise to a negative torque in Fig. [18] agrees very well 
with the isothermal model (see Fig.[T2l. For our nominal density 
(indicated by the dashed line in Fig. [T9l we see that the magni- 
tudes of the torques due to the inner and the outer spiral wave 
are lower than for the low-density case. This is due to compres- 
sional heating in the waves, as can be seen in the top left panel 
of Fig. [14] Because of the higher temperature, the same pressure 
distribution as in the isothermal case requires the density to be 
lower in the waves, and therefore the torque is less strong. 

However, it turns out that the total torque due to the spiral 
waves is still negative. This means that it is the corotation region 
that determines the sign of the torque and therefore the direction 
of migration. In the next section, we will further specify where 
and how this torque arises using local, adiabatic simulations. 

Before we do that, we take a brief look at the accreting 
radiation-hydrodynamical model in Fig. [19] The whole numer- 
ical set-up was the same as for our nominal density model in 
Fig. [19] only now accretion was turned on. All differences be- 
tween the solid line and the dashed line in Fig.[l9]arise from the 
accretion process. Note that the sphere of influence of the ac- 
cretion flow has a radius of approximately 0.1 r p , or 2 H. The 



outflow of energy released by accretion heats the spiral waves 
even more, and consequently the torques decrease in magnitude. 
On top of that, the asymmetric hot bubble seen in Fig. [TBI itself 
exerts a large positive torque onto the planet, which can be iden- 
tified as the large peak in the torque around r = 1 in Fig. [T9j 
This peak is less strong than in the nearly-isothermal case (see 
Fig. [T2l because the heat generated by accretion is diffused by 
radiation over a much larger volume. 

7.4. Local, adiabatic models 

In this section we take a step back and look at local, adiabatic 
simulations. That is, we keep V = 1.4 but we switch off the ra- 
diative transport, and we further reduce the computational effort 
by simulating only a small part of the disk of 5 initial pressure 
scale heights around the planet. This volume is large enough to 
capture the basic features of the flow, because the effects of the 
low-mass planets we consider on the disk remain fairly local. 
We made sure that the global resolution in these local models 
matches the highest local resolution of the previously described 
models. 

In these simplified simulations the gas is still heated by com- 
pression and gravitational effects, while radiative cooling is ab- 
sent. This way, these simulations would apply to a very dense 
region where the cooling time scale is much larger than the dy- 
namical time scale everywhere, also in the upper layers of the 
disk. Although unrealistic for real protoplanetary disks, these 
simulations allow us to pin down where the positive torque in 
a non-accreting model arises while not being hindered by the 
cooling effects of radiation. These cooling effects will partially 
wash out the effects of compression, and while the total torque 
is still positive, its origin is much harder to find. 

In IPaardekoo per & Mellemal (l2006al) . it was found that in 
these models, a warm, underdense plume forms behind the 
planet, which is responsible for the torque reversal. We provide 
some additional discussion here. Further simulations have indi- 
cated that the torque due to temperature effects, 
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is proportional to the radial entropy gradient of the unperturbed 
disk: 
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In our model disk, S is negative, which gives a positive torque 
on the planet. Reversing the entropy gradient, one finds that 
the torque also changes sign. This is also seen using a different 
method (Baruteau & Massetll2007l) . We leave a detailed discus- 
sion of this effect, in particu lar its emergence from linear the - 
ory, to a future publication (Paardekooper & Papaloizou 2007). 
Here, we just note that the mechanism depends on the conserva- 
tion of entropy along streamlines in the horseshoe region, which, 
in the presence of a background radial entropy gradient, leads to 
a strong corotation torque. In a barotropic disk, a similar mech- 
anism operates involving the specific vorticity (see IWarJl 19911) . 
It is interesting to note that, in contrast to the accretion rates, the 
torque displays discontinuous behavior at T = 1, For an adia- 
batic disk, where density and temperature decrease outward, ac- 
cording to Eq. < TT~8T > the entropy gradient reaches a minimum for 
!T — ► 1, leading to a strong positive corotation torque. For V = 1, 
however, we have a locally isothermal disk in which entropy is 
not conserved along streamlines, which kills the entropy-related 
corotation torque. 



8. Discussion 

Our results on gas accretion onto low-mass protoplanets sug- 
gest that mass flow into the Hill sphere depends critically on 
the temperature profile of the planetary envelope. In isothermal 
simulations, pressure gradients correspond to density gradients 
alone, and therefore a pressure-supported isothermal envelope 
will have a very high central density. This leads to the relatively 
high accretion rates in isothermal models. When the energy bal- 
ance is taken into account, part of the pressure gradient that is 
needed to sustain the envelope is carried by a gradient in temper- 
ature and therefore the corresponding density profile will be less 
peaked, and therefore the measured accretion rate will go down. 
Moreover, the energy released in the accretion process raises the 
central temperature even further, slowing down the mass flow 
even more. This way, accretion is a self-limiting process, and the 
subsequent accretion rate is more than an order of magnitude be- 
low the isothermal result. The val ues we find are comp arable to 
the solid accretion rates found by Poll ack et al.l d 1 996b . depend- 
ing on the mass of the planet. Heating due to this planetesimal 
bombardment should lead to a similar drop in the accretion rate 
as observed here. We conclude that accretion rates obtained from 
isothermal simulations strongly overestimate the growth rate of 
embedded planets. 

Because the release of accretion energy plays such a signifi- 
cant role in accr etion, migration and observations of protoplane- 
tary envelopes dKlahr & Kiev 2006|) it is important to consider 
the validity of this accretion recipe. Although at first sight it 
seems a good idea to conserve total energy and momentum dur- 
ing accretion, it is not entirely clear if this is indeed the case. The 
momentum carried by the accreted matter, for example, might be 
transferred to the planet itself, which may induce rotation of the 
protoplanetary core. 

As an example for an alternative to energy-conserving ac- 
cretion, consider an accretion recipe that conserves entropy. 
Because at constant entropy the pressure is related to the den- 
sity as: 



p = K p l 



(51) 



we can, using the ideal gas law, immediately write down that 



T = C p 



r-i 



(52) 



for some constant C. This means that taking away mass from 
the grid actually cools the accretion region. In this case the as- 
sumption is that the accretion energy is transferred to the plan- 
etary core. Note that in using an isothermal equation of state, 
taking away mass from the grid amounts to conservation of en- 
tropy as well. Because accretion plays such a major role in non- 
isothermal planet-disk interaction, obviously in the growth of the 
planet but also in the migration behavior, more work on the de- 
tails of the accretion process is necessary. 

But even for non-a ccreting planet, a hydro static envelope is 
never reached (see also lD'Angelo et al.ll2 003b). Mass constantly 
flows in and out of the Roche lobe, and it is not yet clear what 
role these outflows play in planet-disk interaction. In our simula- 
tions we have found that the outflow is predominantly equatorial, 
irrespective of whether the planet is accreting mass or not. The 
depth of the potential well does influence the outflow, however, 
because it is driven by the high pressure region at the location of 
the planet. We found that this outflow did not alter the torque bal- 
ance in the isothermal limit, but because the gas expands rapidly 
inside the outflow it affects the temperature in the midplane and 
it may therefore be of importance in non-isothermal disks. 



The most interesting result from this study is the observed 
reversal of migration direction with respect to isothermal disks. 
We have tried to isolate the origin of the positive torque arising 
in adiabatic simulations, and we have found that gas heating in 
the region where the horseshoe orbits make their U-turn is the 
dominant effect. This heating is driven by a radial entropy gradi- 
ent, where a negative gradient results in a positive torque on the 
planet. 

Classical corotation torques tend to saturate when the vis- 
cous diffusion time scale ac ross the librating regio n is longer 
than the libration time scale (Ogi lvie & Lubowl2003l) . In our in- 
viscid models, saturation should occur, but a much larger time 
range is needed to see the effect. The computational costs of 
radiation-hydrodynamical simulations are too large to achieve 
this, however. Global adiabatic models may provide more in- 
sight in this matter. The new torques arising from a radial en- 
tropy gradient may be kept unsaturated by energy transport by 
radiation, but only if the disk in radiative equilibrium has indeed 
a radial entropy gradient. More detailed disk models are needed 
to answer this question. 

In our models where we solve the energy equation we have 
focused on a 5 M ffl planet. One run involved a 0.6 M ffi accret- 
ing planet that also moved outward, but more study is necessary 
to establish the trend, in particular whether it is a purely linear 
effect. 

When the planet gets massive enough to substantially change 
the density in its direct environment, the situation changes. Due 
to the lower density near the orbit of the planet the positive 
corotation torques decrease in strength because the cooling ef- 
ficiency increases. Note that the cooling time is proportional to 
the density squared as long as the disk remains optically thick. 
Therefore the outward migration induced by corotational heating 
will quickly lose strength when the planet opens up a gap. For 
gap-opening planets, this mechanism does not operate because 
the disk near the planet can cool efficiently (density reductions 
of two orders of magnitude were measured in 2D simulations). 

This new migration behavior of low-mass planets raises a 
few interesting issues. First of all, because the direction of mi- 
gration depends on the density, there exists a certain radius Rq 
with density p(Ro) where the total torque is zero. This radius is 
the same for all low-mass planets, because it depends on cooling 
properties of the disk only. Inside Rq, cooling is not efficient and 
planets move outward. Outside Rq, planets move inward through 
ordinary Type I migration. Therefore all low-mass planets near 
Ro in the disk will slowly approach Rq. In a simple power-law 
disk there is only one Ro, which, in the Minimum Mass Solar 
Nebula is located around 10-15 AU. When the density structure 
becomes more complicated (for example due to the presence of 
a gap-opening planet) there may be several locations in the disk 
where the total torque is zero. 

As the low-mass planets approach Rq, planet-planet in- 
teractions will eventually come to dominate. When this hap- 
pens depends on the magnitude of the torque in that region. 
Interesting effects w ere observed for dif f erenti ally migrating 
low-mass planets by ICresswell & Nelsonl d2006l) . The planets 
may lock into mean-motion resonances or be scattered out of 
the system. 

As the disk evolves it slowly loses mass. Therefore, in a sim- 
ple power-law disk, Rq will shift inwards slowly. Once a planet 
has reached Rq either by outward or inward migration, it will 
slowly migrate inward due to the dispersal of the disk. Therefore 
the migration time scale for low-mass planets becomes directly 
coupled to the disk life time. While in the isothermal Type I mi- 
gration scenario migration may proceed on any time scale de- 



pending on the disk mass, when the energy budget is taken into 
account this is no longer the case. Low-mass planet migration 
will proceed at time scales comparable to the disk life time, 
which is equal to the viscous time scale of Type II migration. 
Therefore, this scenario of corotational heating provides a solu- 
tion to the problem of fast-migrating cores. 

In using Eq. (Q~|) in calculating the cooling time we have as- 
sumed that the planet resides beyond the snow line. Interesting 
things may happen in regions where the opacity s u ddenl y 
changes, as was pointed out by iMenou & Goodman (2004). 
Also, grain growth affects the opacity to a great extent. A 
disk in which all particles have grown to cm-size will be op- 
tically thin, and therefore the mechanism of corotational heat- 
ing will not apply. However, such efficient growth of parti- 
cles is not consistent with ob servations of protoplanetary disks 
dDullemond & Dominikl I2005T) and therefore some sort of re- 
plenishment of small grains is needed. 

Another important effect that will have to be included is ir- 
radiation from the central star. Although the disk is not severely 
perturbed by the low-mass planets considered in this paper, it 
is desirable to treat the heating of the disk in a self-consistent 
way. However, detailed models including irradiation from the 
central source are computationally very expensive and therefore 
hard to implement in a dynamical model. Flux-limited diffusion 
is not capable of handling possible shadowing effects properly, 
and therefore one has to switch to a ray-tracing method. 

This new migration behavior should in the end b e included 
in pop ul ation synthesis models of planet formation dlda & Linl 
2004a, b, 2005). At the moment, these models require Type I mi- 
gration to be reduced by a factor of a few in order to match the 
observed extrasolar planet population. The new insights in low- 
mass planet migration presented in this paper can increase the 
level of consistency in these models. 

9. Summary and conclusions 

In this paper we have presented the first global radiation- 
hydrodynamical simulations of the interaction between low- 
mass planets and protoplanetary disks. Specifically, we looked 
at accretion and migration rates in the mass regime of linear disk 
response. 

We have found that accretion and migration rates depend 
sensitively on the ability of the disk to radiate away energy gen- 
erated by compression in the various flow regions around the 
planet. Compression deep within the planetary envelope leads to 
a decrease in accretion rate onto the planet of more than an or- 
der of magnitude. Compression within the tidal waves generated 
by the planet leads to a decreased Lindblad torque, which makes 
the planet move inward but at a slower rate than predicted by 
analytical models for isothermal Type I migration. Compression 
in the horseshoe region, finally, leads to a large positive coro- 
tation torque that is able to change the direction of migration, 
depending on the radial entropy gradient. 

All these effects depend critically on the cooling properties 
of the disk. For a density appropriate for the Minimum Mass 
Solar Nebula at 5 AU we find outward migration, while a ten 
times lower density gives rise to inward migration. Lowering the 
density even more we are able to reproduce the analytical Type 
I torque. This means that inward Type I migration is restricted 
to the outer parts of protoplanetary disks, or to older disks that 
have a lower density. 

Because low-mass planet migration is driven by the decrease 
of density in the disk, it proceeds on the viscous time scale, just 
as Type II migration for gap-opening planets. This means that 



corotational heating is able to solve the problem of protoplane- 
tary cores that move inward too fast to form giant planets. All 
planets within the region where the cooling time scale is much 
larger than the dynamical time scale will experience outward mi- 
gration, depending on the local entropy gradient. 

Future studies should be aimed at finding the most appro- 
priate recipe to simulate gas accretion onto the planet, and at 
including irradiation from the central star to calculate the equi- 
librium disk temperature profile self-consistently. 
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